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We study the dissolution of a solid by continuous injection of reactive "acid" particles at a single 
point, with the reactive particles undergoing biased diffusion in the dissolved region. When acid 
encounters the substrate material, both an acid particle and a unit of the material disappear. We 
find that the lengths of the dissolved cavity parallel and perpendicular to the bias grow as i 2 /( d+1 ) 
and t 1 '( d+1 ', respectively, in d-dimensions, while the number of reactive particles within the cavity 
grows as t 2 /( d+1 ). We also obtain the exact density profile of the reactive particles and the relation 
between this profile and the motion of the dissolution boundary. The extension to variable acid 
strength is also discussed. 

PACS numbers: 47.70.-n, 44.35. +c, 05.40. Jc, 64.70.Dv. 
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I. INTRODUCTION 



The dissolution of a solid material by contact with a 
reactive fluid is a fundamental process that underlies cor- 
rosion [Q, diagenesis (|J^1, erosion etching HH, and 
many other industrial processes. The same dynamical 
process can also be viewed as the melting of a solid by 
heating the material at a single point in the interior J7J. 
These types of dissolution (or melting) processes are de- 
scribed by the motion of the interface between the re- 
active fluid and the solid. In many situations, molecu- 
lar diffusion is the transport mechanism for the reactive 
particles, and this leads to diffusion-controlled moving 
boundary value problems MM- 



of the acid particles inside the dissolved cavity, as well as 
the shape and time dependence of the boundary between 
the fluid and unreacted material. 

In Sec. II, we first define the model and write the re- 
action diffusion equation that governs the density of re- 
active particles in the continuum limit. In Sec. Ill, we 
then solve for the steady state density profile of reac- 
tive particles in the dissolved cavity. This profile satis- 
fies the anisotropic Laplace equation, which is the time- 
independent limit of the basic equation of motion. In 
Sec. IV, we investigate the motion of the interface and 
determine the different characteristic lengths of the cav- 
ity in the directions parallel and perpendicular to the 
bias. We briefly summarize in Sec. V and also discuss a 
generalization of the system to variable acid strength. 
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FIG. 1. Schematic illustration of the dissolution process. 
Reactive particles (dots) are continuously injected at rate A 
at a single point (circle). Each particle undergoes biased dif- 
fusion, with bias in the parallel direction. When a particle 
reaches the boundary of the dissolved cavity, a unit of the 
host material and the particle both disappear. 

In this work, we consider the kinetics of this dissolu- 
tion process when there is a superimposed bias on the 
diffusive motion of the acid. Such a bias can be easily re- 
alized, for example, by an electric field acting on ionized 
particles, or by a gravitational field or a pressure gradi- 
ent acting on a flowing fluid. We find that the bias is a 
relevant perturbation with respect to molecular diffusion 
and gives rise to a dissolution process different from that 
caused by isotropic diffusion |i"o| . As a function of time, 
the size of the dissolved region grows continuously and 
preferentially in the direction of the bias (Fig. The 
basic questions that we shall study are the density profile 



II. THE MODEL 

We consider the following microscopic dissolution pro- 
cess (Fig. [l]). Initially all sites are unreacted. Acid parti- 
cles are injected at rate A at the origin of a ei-dimcnsional 
solid substrate. After injection, an acid particle under- 
goes biased diffusion until it hits a site at the interface 
between unreacted substrate and the dissolved cavity. 
In this interaction, both the host substrate site and the 
acid disappear. We can think of the acid as having unit 
strength so that one acid particle and one substrate par- 
ticle are consumed in a reaction. Later we will generalize 
to allow acid to dissolve many substrate sites before be- 
ing neutralized. In the context of melting, we can think 
of particle injection as the localized input of heat and 
dissolution as the melting of the solid when heat reaches 
the interface. 

In the limiting case where the reactive particles un- 
dergo isotropic diffusion, the resulting dissolution pro- 
cess has been extensively studied, both in the context of 
melting (jjand in the framework of diffusion-controlled 
reactions E3|. Here the radius of the dissolved region 
R(t) grows as (tint) 1 / 2 and as t x l d for spatial dimension 
d = 1 and d > 2, respectively. The density profile of the 
acid particles is also radially symmetric and asymptoti- 
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cally approaches a steady state as a function of the scaled 
radial distance r/R(t). 

These results have a simple origin. In d — 1 , for TV dif- 
fusing particles initially located at the origin, the furthest 
particle from the origin after time t will be a distance of 
the order of (tin N) 1 / 2 |0. Thus if At acid particles are 
injected continuously at the origin, the most distant par- 
ticle, and therefore the position of the interface should 
be (tin At) 1 / 2 from the origin. For d > 2, since each acid 
particle dissolves a single substrate particle, the dissolved 
volume can be at most At. Consequently, the radius can 
be no larger than t 1/,d for d > 2. Within the dissolved 
cavity, the density profile of acid particles away from the 
interface approaches a static limit for d > 2. This den- 
sity profile thus obeys the Laplace equation and decays 
as r 2 ~ d . For d < 2 the density profile is not static and it 
can be obtained conveniently by a scaling solution JTTJ]. 

In the presence of the bias, we need to consider sepa- 
rately the growth parallel and perpendicular to the bias. 
The dynamics of this anisotropic dissolution process is 
governed by c(r, t), the concentration of acid at position 
f within the dissolved region at time t. This concentra- 
tion obeys the convection-diffusion equation 
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subject to the absorbing boundary condition c(f, f) = 
for |f| = R(0, t), where R(0, t) is the radius of the moving 
interface as a function of 9 and t (Fig. Q). Here we have 
taken the bias direction as along x. The motion of the 
interface is then governed by the flux of acid onto the 
interface 
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where K is the parameter which quantifies the acid 
strength. Here we define this constant to be 1 and later 
generalize to arbitrary acidity. Notice also that there is 
no convective contribution to this flux (vc) because of the 
absorbing boundary condition. 

A basic feature which simplifies much of the analysis is 
that the density profile of the acid within the dissolved re- 
gion is stationary in time except near the boundary. This 
arises because the input of new particles compensates for 
their loss at the boundary. This same simplifying fea- 
ture, which also applies in the case of isotropic diffusion 
for d > 2, ultimately stems from the transient nature of 
biased diffusion [Q . We now exploit this stationarity to 
obtain the exact concentration profile of acid within the 
dissolved cavity. 



gives the classical Laplace equation with steady-state so- 
lution c ss (r) oc r 2 ~ d for d > 2. To find the corresponding 
solution in the presence of a bias, we perform a Fourier 
transform of the anisotropic Laplace equation to yield 



Dk 2 c(k) + ivk x c(k) + A = 0, 



with solution 
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Inverting this Fourier transform gives the steady-state 
acid concentration 
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In the last step, we complete the square in the denomi- 
nator and then shift k x by k x — iv/2D. The last integral 
in Eq. © is @, 



/ vr \ 
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where is the modified Bessel function. 

This exact solution has very different forms in the 
regions x > and x < 0. In the interesting case of 
a; > 0, we substitute the asymptotic expansion K v {z) ~ 
{■n/2z) 1 / 2 e- z |1 into Eq. @ to obtain 
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In the special cases of d = 1, 2, and 3, this reduces to 
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Conversely, for x < 0, c ss (r) decays exponentially as a 
function of the distance from the origin, with the length 
scale of this decay proportional to D/v. 



III. STEADY-STATE CONCENTRATION 
PROFILE 

Setting the time derivative in Eq. (Q) equal to zero, an 
anisotropic Laplace equation results. For zero bias, this 



IV. INTERFACE MOTION 

To gain a fuller appreciation for time-dependent fea- 
tures and the motion of the interface, we have performed 
Monte Carlo simulations of the dissolution process. Our 
simulations are based on simply tracking the motion of 
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all the reactive particles. Each particle performs a biased 
nearest-neighbor random walk on a d-dimensional hyper- 
cubic lattice, with hopping probability equal to l/(2d) in 
the 2d — 2 directions perpendicular to the bias, and equal 
to 1/ (2d) < p+ < l/d and p- = l/d — p + < p+ in the ±x 
directions, respectively. These hopping probabilities give 
a bias velocity v = p+ — p— , as well as a superimposed 
isotropic diffusion process, with the diffusion coefficient 
in all coordinate directions equal to l/2d. We hence- 
forth fix the injection rate to be A = 1. Each lattice site 
is initially regarded as one unit of solid material which 
disappears when it is contacted by a reactive particle. 

The choice of the bias in our simulations is dictated by 
basic physical considerations. If the bias velocity is too 
small, there is a long crossover time before the bias dom- 
inates over the diffusion. On the other hand, for a bias 
velocity which is close to the maximum value of l/d, the 
length of the dissolved region becomes extremely large 
and this requires considerable computer memory to store 
the data of the system map. For these reasons, we found 
it optimal to consider intermediate values of the velocity. 

As a function of time an elliptically-shaped dissolved 
cavity grow in which the interface remains relatively 
smooth (Fig. ||). Within this cavity, there is a distri- 
bution of mobile reactive particles which have not yet 
reached the interface. These physical characteristics have 
different dependences in spatial dimension d = 1 and in 
higher dimensions; we therefore discuss these two cases 
separately. 
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FIG. 2. Typical shape of the dissolved region on the square 
lattice after t = 10 4 time steps. Reactive particles and sites 
on the dissolution boundary are denoted by crosses and cir- 
cles, respectively. The injection point is at (x, y) = (0, 0) and 
the bias velocity is v — 0.2. For this velocity, vt 3> v2Dt and 
thus the system is far beyond the initial transient regime. 



A. One Dimension 

We first apply a simple flux balance argument jl5| to 
show that for x > the interface boundary R(t) moves 
with a fixed propagation velocity - defined to be vi - 



which is less than the particle velocity v. Since the in- 
put of reactive particles occurs at rate A, the particle 
flux in the +x direction is simply A. Thus a unit flux 
would lead to an interface velocity Uj/A. On the other 
hand, in a reference frame that moves at velocity v, the 
reactive particles are at rest while the substrate parti- 
cles (with density p) move with velocity —v. In this 
moving frame, a flux of substrate particles — pv would 
lead to the interface moving at velocity vj — v. There- 
fore a unit particle flux would give an interface velocity 
(vj — v)/(—pv). Since the reaction has the symmetrical 
stoichiometry A(acid) + B (substrate) — > 0, the two veloc- 
ities under conditions of unit flux must be equal. This 
then leads to the interface velocity 



vi = Xv/ (A + pv). 



(9) 



For r w R, the density profile decreases sharply from 
the constant value X/v (Eq. (||) to 0. Because the disso- 
lution process is equivalent to the reaction A + B — > 
with components approaching each other at finite veloc- 
ity, the width of the reaction front is proportional to D/v 
and does not grow in time |l^,^6[. We have also veri- 
fied these features by Monte Carlo simulation (data not 
shown) . 



B. Dimensions d>2 

For d > 2, let us locate the reactive particles by the 
cf-dimensional cylindrical coordinates r — (x, r±), where 
rj_ is the d— 1-dimensional radial vector perpendicular to 
the x axis. Similarly, we write R — (R\\,R±) to denote 
the position of the interface. Since the dissolved cavity 
grows predominantly along the direction of the bias, we 
focus our attention on this downstream portion of the 
interface. 
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FIG. 3. Plot of Rii versus t on a double logarithmic scale 
for d — 2 and 3 (bias velocity v = 0.3 for both cases). The 
data represent averages over 500 (d — 2) and 1000 (d — 3) re- 
alizations. The inset shows the local slopes of the data versus 
1/lnt. These appear to converge to 2/3 in d — 2 and 1/2 in 
d = 3 (dashed lines), in agreement with Eq. (|l0[). 
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We now determine how Ru and R± depend on time. 
Let us assume that Rn ~ i v n and R± ~ t u± . Since the 
motion of the reactive particles in the transverse direc- 
tion is diffusive, we expect that R± oc R,}/ 2 ~ i^n/ 2 . 
Then the volume of the dissolved region is proportional 
to V ~ - ft^Hi/ 2 . Since V cannot grow faster 

than At, we must have (d + l)un/2 < 1. On the other 
hand, if V were to grow slower than Xt, the number of re- 
active particles in the dissolved region would have to grow 
with time, in contradiction with the steady state density 
profile derived above. Thus V should grow linearly with 
time, from which we conclude that vu — 2/(d+ 1). Hence 

~ < 2/(d+1) R ± ^t 1 ^ d+1 \ (10) 

These predictions are in very good agreement with our 
numerical simulations (Fig. |^) . 
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FIG. 4. Density profile of reactive particles in d = 2 along 
the x axis at equally spaced times on a logarithmic scale. The 
bias velocity is v = 0.2. The data represent averages over 1500 
realizations. Inset: same data on a double logarithmic scale. 
Except for the sharply decreasing interfacial region, the pro- 
file c ss (x,0) oc x~ 1/2 (see Eq. (§)). 

The dependence of iZii on t can also be determined in- 
dependently from the density profile c(r , t) . Similar to 
the case of d — 1, c(r, t) approaches c ss far from the in- 
terface, while c(f, t) assumes a traveling wave form near 
the interface, which rapidly decays from c ss to (Fig. |^). 
From the inset to this figure, we see that the width of 
this front does not grow in time. Using Eq. (^|), we can 
then approximate the equation of motion for the inter- 
face as R ~ c ss /w, where w <~ D/v is the width of the 
front. Substituting the asymptotic expansion for K v (z) 

in Eq. (|J) then gives c ss (i?||) ~ i?jj 1-t ^ . Using this in 

R ~ c ss /w, we obtain Ru ~ (t/w) 2 + , in agreement 
with Eq. @. 

C. Scaling for the density profile 

To obtain the number of reactive particles, it is con- 
venient to express their density profile in a scaled form. 



Based on the time dependences of R\\ and R±, we in- 
troduce the scaled variables £ii = x/t 2 /( d+1 ^ and £j_ = 
r^/t 1 ^ d+1 \ In terms of these scaled coordinates, 

r = (x + r± ) ' 

+ (11) 

Using the asymptotic expansion of K^, and substituting 
the scaled variables into Eq. (||), we obtain the scaling 
form for the density profile 

c(£||,£i,*) ~ 2 exp [-^J^) > 

= H** /(£,[, (12) 

From this form, we easily obtain the time dependence 
of the total number of active particles N(t) to be 

N ( t ) =J d fc(r,t), 

^ t V(d+i). ( 13) 

This prediction is also in excellent agreement with our 
simulations (Fig. ||). 
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FIG. 5. Plot of \nN(t) versus Int. The data are from the 
same simulations as Fig. ^. The inset shows the local slopes 
which appear to converge to 2/3 and 1/2 for d — 2 and 3 
(dashed lines), consistent with Eq. (p"3|). 

Alternatively, N(t) equals the difference between the 
number of injected particles and the volume of the dis- 
solved region. We use this fact to provide a more pre- 
cise form for the time dependence of the dissolved vol- 
ume V(t). By multiplying Eq. (||) by the surface element 
dS, integrating over the dissolution interface, and using 
Eq. (Q), we have 
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d§~ = -D j dS-Vc\ n 
= -D I dV V 2 c 



= — J dV (d t c + vd x c 



X5(r)) 



dN , 



(14) 



The left hand side is simply equal to V. Therefore we ob- 
tain the obvious conservation equation i{V + N) = X. 
Then Eq. (|l|) gives V(t) ~ Xt ~ at 2/(d+1) , where a is a 
constant related to the integral in Eq. (|l3|). 



V. DISCUSSION 

In this paper, we studied the dissolution of a substrate 
when acid particles are continuously injected at a single 
point and there is an external field which causes these 
particles to undergo biased diffusion. The basic quanti- 
ties of interest in this process are the concentration profile 
of the acid and the growth kinetics of the dissolved re- 
gion. Within the dissolved region, the acid concentration 
follows the steady state profile of biased diffusion; this 
is just the solution of the anisotropic Laplace equation. 
The shape of the dissolved region is strongly anisotropic 
with its length growing in time as £|| ~ £ 2 /( d+1 ) while the 
transverse width grows as £j_ ~ 
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FIG. 6. Scaled boundary profiles for acid strength A = 1 
and A — 5 at different times. All graphs are on the same 
scale. The coordinates for A = 5 are divided by (5t) v , with 
v = 2/3 (horizontal scale), 1/3 (vertical scale). The smaller 
contours are for A — 5. 

A simple and relevant extension of our model is to the 
case of variable acid strength. This can be realized by as- 
suming a substrate density p > 1 so that p acid particles 
must hit a given substrate site before it dissolves (weak 
acid) or that an acid particle dissolves A > 1 substrate 



sites before becoming neutralized (strong acid). In the 
current model, p = A = 1. The case A > 1 is equivalent 
to having a particle density in the substrate p equal to 
1/A and with an acid particle dissolving one substrate 
particle. Incorporating this scaling behavior into Eq. (j^) 
for the interface motion we have 



dR 
dt 



D. 
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-z-r — Vc||^ =fl (0 !t ). 



(15) 



For p > 1, the dissolution boundary becomes smoother 
and grows more slowly since it takes many acid particles 
to dissolve each substrate site. For d = 1, this follows 
immediately from the expression Vj = Xv/(X + pv) whic h 
was obtained from the flux balance argument (Sec. IV A| ). 
In general, Eqs. (|l|) and (||) are invariant after normaliz- 
ing c — > c/p and rescaling A — * X/p. This means that a 
change of the substrate density or acid strength will only 
change the time scale through the injection rate. We have 
tested this hypothesis by simulations in which acidity A 
varied between 1 and 160. At short times, the dissolu- 
tion boundary appears to be much rougher for p < 1. 
Asymptotically, it appears that the boundaries for dif- 
ferent values of p approach a common limit. The overall 
effect of varying the acidity is simply to change the time 
scale. However, the subdominant terms to Eqs. ( |io| ) and 
( |l3| ) seem to have strong acidity dependence so that there 
is a long-lived transient correction to this simple scaling 
behavior (Fig. ||). 
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